clc; clear;

%% 参数设置
phi0_list = linspace(-0.8, 0.2, 20);  % phi0 范围
k_list = linspace(0, 0.2, 20);        % k 范围

dt = 0.01;               % 积分步长
T_total = 250;           % 总积分时间
T_trans = 200;           % 暂态丢弃时间
N_steps = round(T_total / dt);
N_trans = round(T_trans / dt);
delta0 = 1e-7;           % 微扰大小

N_phi0 = length(phi0_list);
N_k = length(k_list);

LE_vec = zeros(N_phi0 * N_k, 1);

% 参数矩阵 A, B，固定值
A = [-2.4, 0, 0;
      0, 0, 0;
      0, 0, 0];
B = [1, -4, -3.5;
     0, 1, 2;
    -1, -4, 1.5];

%% 并行扫描 phi0 和 k
parfor idx = 1:(N_phi0 * N_k)
    [i_phi0, i_k] = ind2sub([N_phi0, N_k], idx);

    phi0 = phi0_list(i_phi0);
    k = k_list(i_k);

    % 初始条件
    x0 = [1e-6; 0; 0; 0; 0; 0; 0; 0; phi0];  % 9维初始状态，x110=1e-6，其他0，phi0变化

    % 系统函数句柄，传入k
    f = @(t,x) mCNN(t, x, A, B, k);

    % 计算最大李雅普诺夫指数
    LE_max = estimate_max_LE_Jacobian(f, x0, dt, T_total, T_trans, delta0, N_steps, N_trans);

    LE_vec(idx) = LE_max;
end

LE_map = reshape(LE_vec, [N_phi0, N_k]);

%% 绘制最大李雅普诺夫指数热力图
figure;
imagesc(k_list, phi0_list, LE_map);
set(gca,'YDir','normal');
colormap(jet);
colorbar;
caxis([0 0.15]); % 设置颜色条范围为0到0.15
xlabel('耦合强度 k');
ylabel('初始相位 \phi_0');
title('双参数最大李雅普诺夫指数图 (Fig.6(b1))');

%% 计算最大李雅普诺夫指数函数（基于雅可比扰动轨道法）
function LE_max = estimate_max_LE_Jacobian(f, x0, dt, t_total, t_trans, delta0, N_steps, N_trans)
    h = dt;
    steps = N_steps;
    trans_steps = N_trans;

    x = x0;          % 主轨迹
    dim = length(x0);
    v = delta0 * [1; zeros(dim-1,1)];  % 固定扰动方向

    sum_ln = 0;
    count = 0;

    for i = 1:steps
        % 主轨迹RK4积分
        x = RK4(f, x, h);

        % 计算雅可比矩阵
        J = Jacobian(f, x, 1e-8);

        % 扰动向量线性演化
        v = v + h * (J * v);

        % 每10步归一化并累积log增长率
        if mod(i, 10) == 0
            dist = norm(v);
            if i > trans_steps
                sum_ln = sum_ln + log(dist/delta0);
                count = count + 1;
            end
            v = (v / dist) * delta0;
        end
    end

    LE_max = sum_ln / (count * h * 10);

    % 绝对值大于10设为NaN
    if abs(LE_max) > 10
        LE_max = NaN;
    end
end


%% 雅可比矩阵计算函数
function J = Jacobian(f, x, eps)
    n = length(x);
    fx = f(0, x);
    m = length(fx);
    J = zeros(m, n);
    for i = 1:n
        dx = zeros(n,1);
        dx(i) = eps;
        J(:, i) = (f(0, x + dx) - fx) / eps;
    end
end

